# -*- coding: utf-8 -*-
# Author: Sun Jiawei
# E-mail: sunjiawei@tbea.com
#
"""
ref:
    https://blog.csdn.net/hal_sakai/article/details/51965657/

"""
import pandas as pd
import pywt
import matplotlib.pyplot as plt
import statsmodels.api as sm

from scipy import stats
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.api import qqplot

from visualization import DataSet
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller

"""数据读取
"""
file_path = './data/opsd-time_series-2020-10-06/time_series_15min_singleindex.csv'
ds = DataSet(file_path)
h1, m1 = '1', '0'  # 数据点时刻 1:00
df = ds.get_panel_data('2015', '1', '2016', '1', h1, m1)  # 2015年整年的数据
df['AT_price_day_ahead'].fillna(method='pad', inplace=True)

"""小波变换
使用Daubechies小波,进行一阶分解,获取概貌序列a1和细节序列d1,
"""
wavelet_name = 'db1'
w = pywt.Wavelet(wavelet_name)
a1, d1 = pywt.dwt(df['AT_price_day_ahead'].values, w, mode=pywt.Modes.smooth)
a1 = pywt.waverec([a1, None], w)
d1 = pywt.waverec([None, d1], w)

fig = plt.figure()
ax_main = fig.add_subplot(3, 1, 1)
ax_main.set_title('原始电价序列(2015全年1:00电价)')
ax_main.plot(df['AT_price_day_ahead'].values)
ax_main.set_xlim(0, len(df['AT_price_day_ahead'].values)-1)

ax = fig.add_subplot(3, 1, 2)
ax.plot(a1, 'r')
ax.set_xlim(0, len(a1) - 1)
ax.set_ylabel("概貌序列")

ax = fig.add_subplot(3, 1, 3)
ax.plot(d1, 'g')
ax.set_xlim(0, len(d1) - 1)
ax.set_ylabel("细节序列")
plt.show()

"""差分
从概貌序列的图像可以发现,该序列看上去不够平稳,因此需要将其进行一阶差分
"""
diff_a = pd.Series(a1).diff().iloc[1:]
# 可以对差分后的序列进行adf检验,以判断是否平稳
# adf检验的原假设是序列为单位根,备择假设为序列是平稳的
adf_res = adfuller(diff_a)
print(f"检验的p值为:{adf_res[1]}")  # noqa
if adf_res[1] < 0.05:  # noqa
    print("序列是平稳的")
else:
    print("序列是非平稳的!")

"""ARMA建模
对差分后的细节序列进行ARMA建模
"""
# todo: 模型的p和q的确认仍需要详细研究
#   参考:https://www.statsmodels.org/devel/examples/notebooks/generated/tsa_arma_0.html
plot_acf(diff_a)
plot_pacf(diff_a, method='ywm')
plt.show()

# 建立模型
model = ARIMA(diff_a, order=(3, 0, 0))
# 模型的训练
res = model.fit()
print(res.summary())
# 对残差序列进行检验
# 1. DW检验
residual = res.resid
print(f"DW value:{sm.stats.durbin_watson(residual.values)}")
'''
DW统计量的值越接近2,则越能说明残差序列无自相关性;
DW统计量的值越接近0,则越能说明残差序列正相关;
DW统计量的值越接近4,则越能说明残差序列负相关.
'''
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
ax = residual.plot(ax=ax)
plt.show()

# 2. 是否为正态分布的检验
_, p = stats.normaltest(residual)
if p < 0.05:
    print('残差序列不符合正态分布')
else:
    print('残差序列符合正态分布')

# 3. QQ图
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
fig = qqplot(residual, line="q", ax=ax, fit=True)
plt.show()

# 4. 残差的自相关序列与偏自相关序列图
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(residual.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(residual, lags=40, ax=ax2, method='ywm')
plt.show()

pd.date_range('1/1/2015', '1/10/2016')
